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I. INTRODUCTION 


Real time simulations of aerospace systems have been developed 
at several NASA, centers for the purpose of providing a testbed for 
numerous studies involving use of cocRpit raockups, visual displays, 
pilot/astronauts, and vehicle motion. It is common in real-time 
simulations of large systems to separate the high and low frequency sub- 
systems within the simulation and perform the integrations of the 
subsystems at different calculation rates. This is done to strike a 
balance between accuracy of calculation and capacity of the digital 
computer. Questions arise as to the accuracy of this structure compared 
to single calculation rates and if any interactions arise that cause 
errors that are worse than those expected from an analysis of the sub- 
system above. 

This report describes a study that was done on a linear aircraft 
model that investigates the questions above. Since actual simulations 
are much more complex with many nonlinearities, these results cannot 
be applied directly; however, the study does show where the problems 
are and gives guidelines for selection of sample rates for multiple rate 
simulations. 
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II. BACKGROUND 


The particular system simulations in question typically have a 
fast mode (aircraft short period, dutch roll) and a slow mode (phugoid, 
roll divergence) which can differ by an order of magnitude in the 
respective natural frequencies. These kind of problems are referred to 
as "stiff" differential eq\\ations in numerical analysis circles although 
special techniques to solve stiff equations are not absolutely necessary 
until frequency multiples on the order of 100 or more exist. 

Many different techniques have been reported [l,2,3,4,s] but none 
use different calculation rates for the different modes of the system. 
The reason for the lack of literature (and lack of interest by numerical 
analysts) on multi-calculation rate techniques appears to be due to the 
difficulty in separating systems into fr st and slow modes. Numerical 
integration procedures are never limited to linear systems which are 
really the only ones that can be cleanly separated. 

In aircraft simulations, system descriptions are sufficiently 
close to linear so as to allow separation based on our knowledge of the 
approximate linear version, Furthennore, since the differential equa- 
tions arise from kno\m physical phenomenon which are similar for all 
vehicles and flight conditions, once the separation has been determined 
for one case, it is applied successfully for most others. 

On the other hand, the general problem of analyzing a linear 
discrete system with multiple sample rates has been studied extensively 
[6,7,8,S,io]. Since any numerical integration procedure can be reduced 
to a set of difference equations, and will be linear if the differential 
equations being integrated are linear, these methods are applicable. 
Unfortunately, the methods are very tedious to apply and require large 
amounts of algebra before going to a computer. ■ Application of these 
methods to the aircraft simulation were studied for simple integration 
procedures but judged to be beyond the scope of the study for the more 
realistic and complex integration procedures. 


III. METHOD Of ANALYSIS 


To provide a common yardstick for comparing the various algorithms, 
it was decided to use the frequency response of the aircraft simulations. 
In particular, the transfer function of the longitudinal mode of a DC-8, 
from elevator command to vehicle attitude was selected for study. Two 
methods were employed; 

1) discrete analysis using z-transforms of the single calculation 
rate cases, and 

2) numerical simulation of the multi-calculation rate cases. 


A, The Selected Example 

A DC-8 in approach configuration was selected for study. The 
transfer function between elevator and attitude for this case is [ll]: 


e(s) 

TJS) " 


1.33S 


(s+0.0605)Cs+0. 535) 


(s^+1. 69S+2. S7) (s®+0. 0198S+0. 0267) 


( 1 ) 


which results in the short period and phugoid characteristics as shown 
in Table I. Note the 10:1 difference between the frequencies of the fast 

Table I: EXAMPLE CHARACTERISTICS 



Natural Frequency 

Damping 

Short Period 
Phugoid 

1.62 r/sec (. 258Ha) 
0,164 r/sec (. 0261Ha) 

0.522 
0. 0606 


and slow modes. 

The magnitude and phase of (1) was determined analytically and has 
been included in all the following graphs for comparison (labeled 
"continuous system"). Since this represents the response of an actual 
aircraft with varying frequency of input commands, the goal of all 
digital approximations of this aircraft is to match the continuous 
response as closely as possible, _ 
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Discrete Analysis of the Single Rate Case 

The transfer function in (1) can be written as a set of differential 
equations: 



where for this example: 

— 1,T, 00522 = 0 

a_ = 2.6S13872 y = 6 

^ ia ** 

= 0.09712526 ^ 0 

- 0. 07000953 y^ = W 

K = - 1.338 

= 0. 5955 

bp = 0.0323675 

1 . Euler Vs Integration : 

Euler’s integration [l2] can be simply stated by: 
y(n-hl) = y(n) + T[y(n)] (3) 


Combining (2) and (3) and using first differences to generate Q and 

.. e 

yields: 2 

9(i=) „ ^2 ^+V^^0 

6 (a) " 4 3 2 

e z + 


where 


T 

= 

sample time 

K 

= 

- 1.338 

"o 

= 

bpT^- bj^T + 1 

“l 

- 

b T - 2 

„2 

"'o 


asT - a^T + a^T - 
3 „ 2 

mi 


a T - 2a T ^ 3a T 

A u 

>"2 

= 

a^T^- 3apT + 6 
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The frequency response of this discrete transfer function can be 
determined by evaluating (4) with z taking on values around the unit 
circle. The computer code for doing this is contained in Appendix A and 
the results are contained in the following section for T*s ranging from 
0.05 sec to 0.5 sec. 


2. First Order Adams Integration [is] ; 

The algorithm is; 

y(n+l) = y(n) fef(n) - f(n-l) ] (5) 

A 

where f(n) = yCn) from (2). 

This algorithm makes use of one past value of the derivative func- 
tion and therefore increases the order of the discrete jystem. 

The discrete transfer function of (5) is: 

e Cz) ®4^ ^ ^ ®1^ ®0 

a (z) ^ " 10 ! 9 „ 8 „ 7~~ 5~ 4~ 3 „ 2 „ 7" 

e' + C^z +C„a +C^z +C^z +C_z 4 C.z 4C^Z +C^Z4C„ 

10 9 8 7 6 5 4 3 2 10 

( 6 ) 

where the coefficients are defined in Appendix B along with the computer 
code to evaluate the frequency response for (6). 


3. Second Order Adams Integration [l3]; 


The algorithm is; 

y(n+l) = y(n) + ^ l23f (n)-16f (n-l)+5f (n-2)] 


C7) 


which yields: 
9 (2:) 


14 13 


+ B^z + 


(S> 


'16 


15 


+ C^z + Cq 


where the coefficients are defined in Appendix C. 


05 eoob. 


pA.Glil 
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C, simulation of the Multi-Rate Case 

An alternate way of expressing (1) and (2) is also given ty Toper 

[ll]: 



(9) 


The separation of these equations into fast and slow modes can be done 
in several ways. The ideal manner would be to transform the equations 
into their normal modes which would then produce two coupled 2nd order 
systems, one with pure short period characteristics and one with pure 
phugoid characteristics. Each could be integrated at sample rates 
suitable to that mode. In practice, this is difficult due to the non- 
linear terms in the equations. Furthermore, transforming to and from 
another state definition taltes cpu time and may result in higher cpu 
loading than the more straightforward methods described next. 

1. The 1 X 3 Separation 

If a normal mode analysis of an aircraft is performed, we 
find that the short period consists primarily of or, q, and 0 motion 
with insignificant effect on u. The phugoid consists of u, q, and 
0 with little effect on or. Therefore, since u is the only state 
that does not involve ’'fast" behayior, it is the only state that can 
be safely calculated at the slow calculation rate. The "l X 3" separation, 
recognizes this fact and partitions accordingly. The equations are; 


Fast Loop - 



( 10 ) 
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Slow Loop - 


u ~ - 0.0291U + [0.0629 -32.2] 



( 11 ) 


2. The 2X2 Separation 

Another natural separation is based on fast calculation of 
orientation, q and 8, and slow calculation of translation, u and w. 
It is attractive because a larger portion of the calculations are done 
at a slower rate, hence more cpu time savings appear achievable. The 
equations are as follows. 

Fast Loop - 


• 




*— 

— — 



i— 



■? 

= 

-.792 

0 


q 



-.77 X lO**® 

0087 


U 

bJ 


_ 1 . 

0 

- 

_0_ 


L 0 j " 

0 

0 


w 


( 12 ) 


Slow Loop - 


— ^ — 





— ~ 


— *— 


—* 

— 


"• — 

u 

= 

-.0291 

.0629 


u 


0 


0 

-32.2 


q 

w 


_-.251 

-.628 _ 


w 


_-10.2_ 

e 

_243.5 

0 


_e_ 


(13) 

3. Simulation Procedures 

The frequency response of each separation was determined using 
Euler's Integration (3) and the 1st order Adams Integration (5). It 
was evaluated at calculation rate ratios (IR) varying between 1 and 20, 
Since IR - 1 is the single rate case, these calculations could be 
checked by comparing with the analytical evaluations described in III-B, 
The frequency response was determined by evaluating equations 
(10) and (11) or (12) and (13) using the integration formulas with 

equal to a sine wave of magnitude = 1. After an initial transient 
settling delay the resulting sine wave magnitude and phase was assumed 
to be the desired frequency response. The short period portion of the 
transient response was quite short; however, the phugoid transient 
response was unduly long to wait for settling. Therefore, the procedure 
calculated the phugoid transient response based on the continuous system 
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(1) using inverse Laplace transforms and subtracted this from the 
numerical evaluations before determining amplitude and phase# Appendix 
D contains a listing of the computer code. 


IV. RESULTS 


Table II contains a summary of the figures which represent the results 
of this study. In the simulations, all IR's between 1 and 20 were eval- 
uated. Those cases not shown in the figures were found to be unstable. 


TABLE II: SUMMARY OF FIGURES 


Figure 

Analytical 

Simulated 

Integration 

Algorithm 

T 

fast 

. IR 

Separation 

mm 

X 


E 


1 

— 

B 

X 


A1 


1 

— 


X 


A2 


1 

— 

B 

X 


E,A1,A2 

.05 

1 

— 

5 

X 


E,A1,A2 

.1 

1 

— 

6 

X 


E,A1,A2 

.2 

1 

— 



X 

E 

.05 

1-20 

1X3 



X 

E 

.1 

1-10 

1X3 



X 

E 

.2 

1-5 

1X3 

10 


X 

E 

.05 

1-10 

2X2 

11 


X 

E 

.1 

1-5 

2X2 

12 


X 

E 

.2 

1-3 

2X2 

13 


X 

A1 

.05 

1-20 

1X3 

14 


X 

A1 

.1 

1-10 

1X3 

15 


X 

A1 

.1 

1-5 

1X3 

16 


X 

A1 

.05 

1-10 

2X2 

17 


X 

A1 

.1 

1-5 

2X2 

18 


X 

A1 

.2 

1-3 

2X2 


The most significant result is the difference between the two 
separations. This can be seen by comparing the deviations from the 
continuous curves in Figs. 7, 8 and 9 with those in 10, 11 and 12 
respectively for Euler's Integration and similarly Figs. 13, 14 and 15 
with 16, 17 and 18. For both integration methods, the 1 X: 3 separation 
is decidedly superior. This is no doubt due to the fact that the 2x2 
separation solves the w< or a) equation at the slow rate while this 
state is important to the short period dynamics. 
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The first order Adams integration appears to be the best choice 
of Integration methods. Its advantage over Euler's method and small 
disadvantage compared to a second order Adams is best Illustrated in the 
Fig. 6b phase plot; however, examination of the magnitude in Fig. 6a 
shows the Euler method's error arising at a lower frequency but remaining 
smaller at the higher frequencies. The same kind of behavior is exhibited 
at the faster sample rates (Figs. 4 and 5 ) but is more difficult to see 
due to the greater accuracy. 

The sample rate requirement for aircraft simulation with a first 
order Adams integration is dependent on the desired input frequency to 
be adequately simulated. Examination of Figs. 13, 14 and 15 indicate 
that one should select the fast sample rate at approximately 10 times 
the input frequency to be followed and that a slow rate at one-tenth 
this rate yields no degradation. In other words, to follow a 2 Hz input, 
one should solve the short period equations at 20 Hz and the phugoid at 



V. CONCLUSIONS 


For a linear model of longitudinal aircraft motion, separation of 
the equations of motion into slow and fast calculation, rate groups is 
best accomplished by performing u integration at the slow rate and 
w,q,0 at the fast rate. A separation with u and w as the slow 
variables and q and 0 as the fast gave substantially less accuracy. 

A first order Adams integration procedure appeared to be a good 
choice for real time aircraft simulations. 

For the example used (w , ^ , , *= 0. 25 Hz, oj , , , ~ 0.025 Hz), 

short period phugoid 

the fast sample rate should be selected at approximately 10 times the 
maximum input frequency for which accurate aircraft simulated response 
is desired. A slow rate of one-tenth the fast rate yielded no degrada- 
tion over the single rate case. 
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APPENDIX A 


FREQUENCY RESPONSE EVALIATION OP EQUATION (4) 
(EUIiER'S INTEGRATION) 


The follovdng computer cod© was used to analytically evaluate 
a G(z) by letting z take on values around the unit circle. The 
coefficients stored Al, A2,... and B1»B2 are for the continuous 
transfer function , equation (1) that represents the example described 
in this report. Lines 21 through 29 compute the coefficients required 
for an Euler Integration of the continuous system. ' 


c cizireM A I'KAhsi'sa ruNcxiSii isi tai A-eUEie cu tu£ tsunns^ fitii 


U ll-l 

C B(ti) i » D{S-1) i. ♦ i f » d(1) 
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* CO 


C flits iBCCHAa CAi.-wULA7££i, tAuU1.’.V£? A.41 li't.’nS £sie /LK^U£ICT aC:>e3>'Se 
C Qk' titli SYn£B ilhCSI! raA^3?ES fUitCX'IZi 19 UlVSti ti( tUC rOBS, 
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... — 6 .. 
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3_ 

10 

11 

^12. 
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11 32. .ilJ jK1jX2. 

i,9i>XCAl. i.i;532 (21) . fa33Eti<1ZI) 


017 A Al.AZ,Ai.Al/n 7 1 )} 622 b), 2 .ji&l 3 «> 72 {»O,C. 0 !> 7 t 3326 f>llt>. 0733 69590 /. 
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APPENDIX B 


FREQUENCY RESPONSE EVALUATION OF EQUATION (6) 

(1st ORDER ADAMS) 

The computer code in Appendix A was modified by replacing lines 21 
through 29 with the following coefficient evaluations. 


21 

CONS 

= -1.336 

22 

Cl * 


23 

^ fO - 

-1Gtt*-i*tiU 

24 

C3 » 

54*'i'#¥4 

25 

c4 = 

-12»r*^4 


C5 = 

27 

Cl - 

16. 

26 

D?. = 

24*fl1*I-64 

23 

C3 = 


30 

04 = 

S4*.\3*rv»3-&4*&2*'i»*2*36<»Al*I-a4 

31 

DB = 


32 

J1S_=- 

-rJi) tt* A4 f.£.t t3^Lj2*i2* i **2 ►.8A&Ltl 

33 

D7 = 

54*A4''I*<‘4-20^Al»l*‘’‘‘3*4*A2*t»*'2 

34 

D0 ® 

-12*ftfe*r**4>Z*A3*2**3 
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ng = 



36 

til = 

4. 

37 

(i2 >= 

6*6l>i'-9 


t; ^ a 


39 

«4 a 


40 

05 = 

62»‘X‘'*2 

ai 

Ki j 

, . . ..... . ... 

42 

A2 ^ 


43 

K3 = 


4i( L- 

L0_=f- 

C'V*i:*i -- - - 

45 

B(1) 

c4*Oi5+CB*g4 

46 

a (2) 

* C3»o5TC4*G4»CS'‘tJJ 

47 

H / 31 


46 

a (4) 

= Cl*&StC2*‘G4*C3*03tGfH'G2<-C5*31 

49 

6(5) 

= Cl*GtM-C2*G3*C3*02*C!l*Gl 


fi (4 V 

= l*i.2+c.2*:: 2i-r! 

SI 

■ 

B(7) 

= Cl*'G2*C2*‘Gl 

52 

6(6) 

c el*Gl ‘ ' 

4? 

CO— a. 

nguK’'^ 

54 

C(l) 

S 08*K3*D3*K2 

55 

C(2) 

? D7*‘K3*D6*it2*D9>Al 



f<3V ^ fif.*k34.m*i{2»nfl*k1 

57 

C(4) 

* D5*K3+DS*K2tli7*'Kl 

56 

C(5) 

* C4»K3*C5*R2m6*it1 

—55 

C4&.)- 

= D3.*K3tU4.+K2H)Bti{1 

60 

C(7) 

= D2♦^3+b3^i^2^04^K1 

61 

C(6) 

s Cl*K3t02^K2+03»>M 



63 

C(10] 

1 * p1«K1 ' . 


52 



APPENDIX C 


FREQUENCY RESPONSE EVALUATION OF EQUATION (8) 
(2nd ORDER ADAMS) 


The computer code of Appendix A was modified by replacing lines 21 
through 29 with the following coefficient evaluations. 
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C3 a 
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54 
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57 
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S3 
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il (5) c 

62 

B(6)x 

63 

E(7) = 
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73 

. . -C(2) = 
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77 
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78 
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C(13) 

65 
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APPENDIX D 


tion. 

1X3 


SIMUIATION FOR FREQUENCY RESPONSE EVAEIATIONS 
IN THE MULTI-RATE CASES 

The following code performs the calculations using Euler’s Integra- 
Note that lines 1G7 through 112 are shown twice, once for the 
separation and once for the 2X2 separation. 
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The modifications to the previous code so as to use first order 


Adams integration are shown below. 
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